Self Diffusion in Simple Models: 
Systems with Long Range Jumps. 

by 

A. Asselah 
Department of Mathematics 

Rutgers University 
New Brunswick, NJ 08903 

R. Brito 
Departamento Fisica Aplicada I 
Facultad de CC. Fisicas 
Universidad Complutense 
28040-Madrid, Spain 

and 

J.L. Lebowitz 
Department of Mathematics and Physics 
Rutgers University 
New Brunswick, NJ 08903 



2 

We review some exact results for the motion of a tagged particle in simple models. 
Then, we study the density dependence of the self diffusion coefficient, Djv(p), in 
lattice systems with simple symmetric exclusion in which the particles can jump, 
with equal rates, to a set of N neighboring sites. We obtain positive upper and 
lower bounds on F N (p) = N((l - p) - [D N (p)/D N (0)])/(p(l - p)) for p E [0, 1]. 
Computer simulations for the square, triangular and one dimensional lattice suggest 
that Fn becomes effectively independent of N for N > 20. 

1. Introduction 

Many properties of macroscopic systems are universal, retaining their qualita- 
tive features under drastic simplifications of the underlying microscopic structures. 
Thus, lattice gas models have greatly enhanced our understanding of phase tran- 
sition phenomena in equilibrium systems. Their dynamical behavior, currently an 
active area of research, promises to be similarly fruitful for understanding nonequi- 
librium properties of macroscopic systems. 

This article explores some aspects of self-diffusion in lattice models. After a brief 
overview of some rigorous results, we derive new results for systems with long range 
jumps. It is dedicated to Matthieu H. Ernst, a leader in the field of kinetic theory 
and lattice gases, on the occasion of his sixtieth birthday. 
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interacting particle system. A tagged particle is exactly like any other particle in 
the system, its 'tag' permits us to follow its trajectory X{t). This yields a relatively 
simple probe of time correlations in a system of interacting particles in an overall 
stationary state. 

The self-diffusion coefficient D s is defined, in an infinite stationary system with- 
out drift, as [1] 



where d is the spatial dimension of the system and the average () is over the 
stationary measure. We expect that in a real fluid the limit (1.1) will exist, be 
positive, and be given by the Einstein-Green-Kubo formula 



where (v(r)-v (0)) is the velocity autocorrelation function [1]: A simple computation 
gives ((X(t) - X(0)) 2 ) = 2 f*(t - t){v(t) ■ v{0))dr, so (1.2) reduces to (1.1) when 
(v(t) ■ v(0)) decays sufficiently rapidly. 

The self-diffusion coefficient is a global dynamical parameter associated with 
macroscopic system in equilibrium, i.e. spacially uniform. Therefore, it is gener- 
ally different from the bulk diffusion coefficient Dj,, which relates to the evolution 
of a nonuniform density in a non-stationary system. D s can be thought of as a 



D s = ±- lim h{X(t)-X(0)] 2 ), 
Za t^oo t 



(1.1) 




(1.2) 
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components of a system which differ only by a property, say color, that plays no 
role in the dynamics, while the overall system, ignoring color, is in a uniform state 
[1, 2]. An approximate experimental realization of such a situation occurs when 
the components are isotopes of He 3 atoms whose spins are polarized in different 
directions. 

Going beyond (1.1) and (1.2), we can introduce a "scaling" parameter e and 
define X e (t) as e[X(t/e 2 ; •) — X(0)], where the • indicates the dependence of the 
trajectory on the coordinates and velocities of all the particles at t = [1,2]. 
Typically, we expect that in the long time limit, e — > 0, the process {X e (t), t ElR} 
converges in probability -after subtracting out any drift- to the law of a Brownian 
motion {Wo s (t), t G M} with diffusion coefficient D s given by (1.1) [1,2]. We 
summarize this by 

UmX e (t) = W Ds (t). (1.3) 

The behavior (1.1)-(1.3) has been proven for the one component, one dimensional 
system of hard rods with diameter a [1, 3, 4]. For this idealized system D s can be 
computed exactly, 

D s {p) = ^—^-{\v\), and (H)= f \v\h(v)dv. (1.4) 

P J-oo 

Here h(v) = h(—v) is the one particle velocity distribution function; this need not 



that (1/p — a) is the mean free path in this system, the interpretation of (1.4) is 
very simple. On the other hand, the velocity autocorrelation function (v(r)v(O)) 
depends non-trivially on h(v): it decays like an exponential when h(v) vanishes 
near v = 0, and like t~ 3 when h(v) is Maxwellian [5]. 

The only other continuum system for which the existence of the limits (1.1)-(1.3) 
has been proven is a system of interacting Brownian particles ([6,7]) which models 
suspensions of polymers or even of small macroscopic balls in a fluid. Actually, one 
needs to assume ergodicity of the dynamics, and formula (1.2) has to be modified 
because instantaneous velocities are no longer well defined [1]. 

In d = 1, stochastic models in which the particles cannot cross each other behave 
differently from the mechanical model which yields (1.4): (X 2 (t)} ~ \ft so D s = 
[1,8]. Interestingly, however, yftX (t / 'e 2 ! ; •) still goes to a Gaussian process; see [4] 
for a simple derivation of the one dimensional results. 
2. Lattice Models 

a) General Dynamics 

We consider now systems with one type of particles whose total number is the 
only quantity conserved by the dynamics. We expect however that much of our 
discussion will remain valid also for models where momentum is also conserved 
[7, 9]. The stochastic dynamics of these systems consists of particles "jumping" 



a rate c(x, y; rj), where rj is the configuration of the system just prior to the jump: 
rj = {rj(z)}, with r)(z) = 0, 1, 2, .., specifying the number of particles at site z. We 
shall generally consider the <i-dimensional (simple) cubic lattice Z d . 
The system will be in a stationary state with measure v whenever 

c(x, y; Tj)v(rj) = ^ c(x, y; rj x > y )v(rj x > y ), (2.1) 

x, y x,y 

where rj x,y is the configuration which arises from rj after a particle has jumped from 
x toy. A simple way to satisfy (2.1) is to have the equality hold for each term in the 
sum. The rates are then said to satisfy detailed balance with respect to v. In such 
cases v can be written in the form of a Gibbs measure, v eq {rj) ~ exp[— j3H{rj)], where 
H(rj) is the energy of a configuration rj, and (3 is the reciprocal temperature: see 
[10] for a detailed discussion of Gibbs measures. Detailed balance then corresponds 
to 

c(x, y- V )/c(y, x; r,^) = exp{-/?[#(r/^) - H( V )}}. (2.2) 

In the probability literature a stationary process whose rates satisfy detailed balance 
is called reversible: a film of the system in the stationary state would look the same 
if run backwards. 

The trajectory X(t) now takes values on the lattice. However, after scaling with 
e and letting e — > 0, as in (1.3), the limit will again be a continuous process. 



One of the simplest dynamics for a system of particles on a lattice is the so called 
"zero-range" process [11]. This corresponds to the jump rates c(x,y;r}) depending 
only on rj(x), the number of particles at site x, 

c(x,y;r]) = \g(ri(x))p(y-x). (2.3) 

Here A is a constant, A > 0, while g and p satisfy the conditions 

<7(0) = 0, g(k)>0, for k>0, (2.4) 



p(0) = 0, p(r) > 0, P( r ) = L ( 2 - 5 ) 

rei d 



The stationary measures for this dynamics in the macroscopic (infinite volume) 
limit, are a translation invariant family of product measures, v p , parametrized by 
the average density p[ll]. The probability of having exactly j particles at any given 
site is 

Wj = W) W(U j = Q > 1 > 2 >- ( 2 - 6 ) 

where 

G(0) = 1, G(j)=Ul ig (l), j>1, (2.7) 

and the parameters b and Wq are determined by the normalization and the specified 
average density, p > 0, i.e. 

oo oo 

^W, = l, and ^2jWj=p. (2.8) 



An easy check shows that these measures satisfy the detailed balance condition 
(2.2), with (3H{rj) a sum of single site energies equal to — logWj, if and only if 
p(r) = p(—r). 

Two particular cases of the zero range process deserve mention. When g(l) = I 
the dynamics corresponds to that of independent particles. This gives rise to the 
Poisson distribution 

W S = (p J /jl)e~ p , (2.9) 

Taking g(l) = 1 — <5o,z, corresponding to only the 'top' particle jumping, yields a 
geometric stationary distribution 

Wj = -^—(-^—y. (2.io) 

The stationary measure seen from the tagged particle is the 'Palm measure' 



7/(0). 



-V 



p- 



As the "waiting time" of the tagged particle depends on the number of particles at 
the same site, its average jump rate is given by 



A = V-^ W k , with X k = Xg(k). (2.11) 



Let the displacement, after K steps, of the random walk specified by transition 
probability p(r) be X K . Assuming for simplicity that there is no drift, 

TTt( t\ — n ( 9 19^ 
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we have 

(X 2 K+1 ) = ((X K + Y K f) = {Xl) + (Y£), (2.13) 
where Y k is the displacement of the particle at the (K + 1) step. Clearly, 

(>£) = 5> 2 p(r) = 2D , (2.14) 

and 

(X 2 K ) = 2D K. (2.15) 

A little thought shows that for the zero range process, the trajectory of the tag 
will look the same as the trajectory of a single particle performing a random walk 
on the lattice with transition probabilities p(r). The only difference is that the 
"waiting time" at any site will generally depend on the number of particles there. 
In fact [9], As soon as the process is ergodic, the scaled trajectory X e (t) satisfies 
(1.3) with 

D{p) = \(p)D . 

Ergodicity is easy to show for g(k) = k, and was shown in the case g(k) = 1 — Sk,o 
in [12]. In fact, (1.3) is proven in [13] for all g(k). 

For the case of independent particles, A is independent of p and equal to A, while 
for g(k) = 1 - 5 fc ,o, 
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so that D(p) decreases with density for this case. The opposite behavior is clearly 
also possible. 

Looking back on our arguments leading to (2.15) we see that the main ingredients 
are the independence of the step Y K from the past history of the process (e.g. X k 
is a martingale). This means that (2.15) should remain valid whenever 

c(x,y;ri) = h(jr, x)p(y - x) , (2.16) 

i.e. as long as c(x, x + r; rj)/c(x, x + r'\ rj) is independent of r\ (and by translation 
invariance of x). 

A particular example of (2.16) is (a generalization of) a model due to van Beijeren 
[14], in which 

k 

c(x, y;r]) = Y^ *i9i(jl(x))pi(y - x), 
i=i 

with the <7j and pi satisfying the conditions (2.4), (2.5) and (2.12). The stationary 
measure is now not known in general. In fact we expect that it will have very long 
range correlations [14], yet (1.1) and (1.3) should still be valid with 

i 

where A(p) is the average rate in the stationary measure. 

Remark: It is clear that the diffusion constants D and D are, for anisotropic 
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When (2.12) holds and p(r) is isotropic with respect to the lattice directions then 
D is diagonal. 

c) Models with Exclusion 

We consider now the case where there is a hard core interaction between the 
particles, forbidding the presence of more than one particle at any lattice site. The 
configurations of the system are now given by rj = {r](x)} with rj(x) = or 1 and 
c(x, y; rj) = when rj(y) = 1. The simplest dynamics for these systems correspond 
to the so called simple exclusion processes, in which the jump rate from a site x to 
a site y is independent of the configuration at other sites of the lattice, 

c(x, y; rj) = At/(x)(1 - r]{y))p{y - x), 

with p(r) satisfying (2.5). 

The translation invariant stationary measures v p , with < p < 1, will again be 
product measure with 

W = l-p, W 1 = p, Wj = for j>2. 

The Up will, like before, satisfy detailed balance if and only if the jump rates are 
symmetric, p(r) =p(—r). 

The behavior of a tag trajectory X{t) is now considerably more complicated 
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will influence the probabilities that certain sites are empty and hence the future 
position of the trajectory. 

The proof that (1.1)-(1.3) hold for these models was given by Kipnis and Varad- 
han [5] for the reversible case: if d > 2, then D(p) > when p < 1; if d = 1, then 
D{p) — unless p(r) > for \r\ > 1. For the non-reversible case satisfying (2.12) 
the result is due to Varadhan [15]. As already mentioned, in d = 1, with jumps 
limited to nearest neighbor sites, p(l) = p(— 1) = |, the mean square displacement 
only grows like Vi [8]. 

In the asymmetric case, with p(l) = p, p(— 1) = 1 — p, p > |, (1.1) still holds 
after subtraction of the drift, i.e. for X{t) = X(t) — vt , where the velocity v is 
given by, v = (2p — 1)(1 — p). This was proven for d = 1 by Kipnis [16] and for 
d > 3 by Varadhan and Yau [17] who also prove (1.3) (there is no proof for d = 2). 
Somewhat surprisingly the diffusion constant for X(t) in d = 1 is equal to the drift 
£>(p) = (l-p)(2p-l) [16]. 

The dependence of D on p is not known even for the simple exclusion process. It 
is intuitively clear that D(p) — > as p — > 1. Varadhan [18] proved that the so-called 
"correlation factor" D(p)/[1 — p] is a decreasing function of p bounded away from 
zero as p — * 1. This confirms the behavior found in numerical results for nearest 
neighbor jumps [19]. He also showed that D(p) is continuous in all dimensions and 

j-U„-i- e — „"\ :„ t : u:j-„ , I -nr „\ n/ „/M ^„l„ „/l 
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Remark: We remark here for completeness that many authors (beginning with 
Einstein) have studied situations where there is a special particle with a different 
dynamics than the other particles of the system: in particular, the case where an 
external field acts only on the special particle (see [1]). Though we will not discuss 
this problem here in any detail, we want to mention that recently, Landim, Olla 
and Volchan [20] have studied a one dimensional system where the special particle 
jumps with probability p > 1/2 to the right and 1 — p to the left -with an exclusion 
rule- while all the other particles follow a symmetric dynamics with p = 1/2 (recall 
that the self diffusion constant is in this system [8]). They showed that X e {t) 
converges in probability to a number v(t) which solves a differential equation and 
depends on the initial macroscopic density profile. Their result holds for a large 
class of initial profiles. For instance, when the initial measure is a product Bernouilli 
measure of density p, they showed that 



3. Long Range Jumps 

We discuss now the situation where particles make jumps to a symmetric neigh- 
borhood U, containing N sites, with equal probability, p(r) = iV _1 for r e U , 
p(r) = otherwise. We shall be interested in the behavior of the diffusion constant 

7-1 f „\ „„j a„„.,:i„u..„„„i„j (-„„: — + — vAT/jA „.i at i ™„„ i„ ,„ t„+. + 




l-p 



P 
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as N increases the tag is less likely to revisit, during a fixed number of jumps, a 
previously occupied site and hence there will be less and less memory left of previ- 
ous jumps. The only effect of the hard core exclusion will then be to slow down the 
jump rates by a factor (1 — p), the fraction of attempted jumps which are unsuc- 
cessful due to the target site being occupied. Intuitively this will lead to a density 
independent correlation factor, D N (p) /(D N (0)(1 — p)), in the limit N — > oo: this 
limit is analogous to the van der Waals or mean field limit in equilibrium systems 
when the particles interact via a long range Kac potential [21]. 

Since our dynamics is reversible the result of Kipnis and Varadhan [6] applies, 
so that for each N 

lim^= = WM^ (3-1) 
^D N (p) 



and one expects that 



lim ^44 = 1 - P- (3-2) 

N^oo D N (0) ' V ' 



Actually, we show more: 

Proposition. Set Cjv(p) = Dn(0)(1 — p) — Djq{p)- Then there are positive 
constants c\ and such that 
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In order to determine the behaviour of -Djv(p), we have performed numerical sim- 
ulations in one and two dimensions with periodic boundary conditions at different 
densities. The average of [X N (t)] 2 over many realizations (from 100 to 1000) was 
plotted against t and fit to a straight line passing through the origin. The fitted 
slope is then taken for the diffusion coefficient. 

The one dimensional lattice had 2000 sites. Typically, 200 realizations were 
run for 1000 time steps. The maximum N considered was 100. In two dimensions, 
square and triangular lattices were used, with 200x200 lattice sites. 300 realizations 
were run for about 500 time steps and N < 90. Simulations with larger N were also 
run, but the statistical accuracy was not enough to extract any information beyond 
the zero order one. In Fig. 1 we plot Fn(p) vs. p for the square, triangular, and one 
dimensional lattices for different values of N: the set U to which the particle jumps 
being the N closest neighbors as measured by the number of bonds required to 
reach the site. It appears that F^{p) varies linearly in p with a slope independent 
of dimensions. 
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Appendix. Proof of the Proposition 

We will work in the moving frame of the tag particle and for simplicity write the 
proof for cubic lattices. Thus, U = {y ^ 0: \ \y\\ < n} and N = (2n + l) d — 1. The 
generator of this process is —A with 



where T xy r\ = rf ,v , f y shifts the configuration by a vector y, and we have denoted 
fy(T 0y ) by Ty. The process seen from the tag particle is reversible with respect 
to Opij]) = v p (r]\r](0) = 1) with p in [0,1]. The expectation with respect to P p is 
denoted by E p . 

It is well known, see [1] , that if St is the semigroup generated by A 



-Dtv(o) = — Yl yl~ n2 and w iv) = T7 Yl yi(v(y)-v(y)), 



where y = (yi, . . . , yd), y = (— yi, y2, ■ ■ ■ , yd) and because the diffusion matrix is 
diagonal, we chose w to be the current in the e\ direction. Note that Cn(p) = 



yeu 



a;€^ d \{0} yEx+U 




\\v\\<n,yi>0 



f (Stw, w)dt. 



We introduce the normalized variables 



(v(y) - p) 



such that 



E p [r x r y ] = 5 : 



\Jx,yEZ d \{0}. 
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Recalling a variational formula used in [17] 



we just need to show that for any local function /, 



EP{wf] 2 <cp{l-p)^- d EP[fAf]. 



Now, 



\\y\\<n,yi>0 



Y, V^{r v (T- v , yS - })]. 

||j/||<n,yi>0 



First, we fix (y^, ■ ■ ■ , y°) and work on the line y = (yi, y%, ■ ■ ■ , y°) for y\ G [—n, n]. 
For each y\ > 1, there are y\ — 1 different ways of joining y and y in three steps 
(y, kei + y, (k, y%, ■ ■ ■ , y) with = 1, . . . , y\ — 1 while remaining on the line 
joining y and y. It is important to note that kei+y is at y\ units from (k, y%, ■ ■ ■ , 
and thus as y\ ranges from 2 to n, a given pair of sites on the same line will be 
used at most once. Now, with the notations k\ = ke\ + y and k^ = (k, y%, ■ ■ ■ , y^) 

Ty,y ~ I = Ty ) k 1 Tk 1 ,k 2 Tk2,yTk 1 ,k 2 Ty,k 1 ~ Tq 



= T1T2T3T4T5 - T = t T o • ■■T m -i](T m - T ). 
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Also, the product measure being invariant under exchanges 

Vm G [1,5] EfryiToT, . . .T m _x](T m / - /) = E^r ym (T m f - /), 

where y m belongs to {y, k u k 2 , y}. 

Now, if yi, y 2 belong to a line parallel to ei, we want to 'split' T yijV2 -which arose 
in our previous decomposition- into T yi;Z and T ZjV2 where z is a common neighbor. 
yi and 1/2 will have of the order of n d common neighbors if ] | y± — y 2 | | < n. However, 
for each neighbor, say z, the pair (y l5 z) will be used at most 2n times because in 
our line, y\ has 2n neighbors at most. Thus, we end up with 

J2 yiEO\r y (Ty, y f - f)\ < J2 EO\r,(T x , y f - f)\ 

\\y\\<n, yi >0 ' \\x-y\\<n 

\\y\\<n 

where a depends just on the dimension d, and 7 is an inocuous index different for 
each pair (x, y). 

Denoting by T m a generic exchange operator and 7 an arbitrary index and using 
the Cauchy inequality 

1 t?p r~. <<t\ r j?m 1 ^ / T7ipr„2i T^p\(rr\ r ^\2i\l/2_ fW^TFr 1 rvx\ 
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gives 



\\x-y\\<n 
\\v\\<n 



\x-y\<n 
x,y^O 



2 

n 

<c-p(l-p)E?[fAf}, 



and the first inequality follows. 
b)C N (p)>cp(l-p)g. 

Choose / = E|| y ||<n, yi >o r ^ and write 

0<\\y\\<n \\x-y\\<n 

x,y=lO 

The reason to choose such an / is that 

2/1=1 

Now, r x r y is equal to r x+2/ if y ^ —x and to r x if y = —x, so, for any y, it is easy 
to see that 

E%Tyf-ff]<N. 

Thus, 



E £ P [(rJ-/) 2 ]<iV 2 . 
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Now, it is also easy to see that 



\x-y\\<n 



EP[(T x , y f-f) 2 ]<N 2 , 



and thus 



which completes the proof. □ 
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Fig.l Simulation results of the quantity F N (p) versus the density p. Results of 
simulations in one dimensional lattices are represented by circles, in a square lattice 
by squares and in the triangular lattices by triangles. The number of neighbours is 
shown in the legend. The slope of Fpj(p) as a function of p appears to be independent 
of the dimension or the type of lattice. The intercept, on the contrary, appears to 
depend on the dimension but not on the lattice structure. 





i 1 * ■ 



1 


2 S 




• N=20 


• 


♦ N=30 




N=40 


■ 


■ N=20 




■ N=28 




^ N=80 


A 


a N=30 


A 


a N=54 


▲ 


^ N=96 



0.0 



0.2 



0.4 



0.6 



0.8 



1.( 



